/**
 Fint.h

 \file
 Declaration of class Fint.

 \author Thomas Burton 
 \date 3/13/12
 \copyright 2012 BNL. All rights reserved.
*/

#ifndef _FINT_CERN_E104_CPP_
#define _FINT_CERN_E104_CPP_

namespace tppmc {
   
   /////////////////////////////////////////////////////////////////////////////
   // C++ implementation of the CERN E104 routine FINT for multidimensional
   // linear interpolation.
   // http://dollywood.itp.tuwien.ac.at/cernlib/shortwrupsdir/e104/top.html
   // Arguments follow the same order as the Fortran FINT subroutine.
   /////////////////////////////////////////////////////////////////////////////
   
   template<typename T>
   class Fint {
      
   public:
      
      // Note that the last argument differs in format from the Fortan version.
      // Instead of a multidimensional array, the array is one-dimensional.
      virtual T operator()( const int nDimensions,
                           const std::vector<T>& interpolationPoint,
                           const std::vector<int>& arrayDimensions,
                           const std::vector<T>& knownPoints,
                           const std::vector<T>& valuesAtKnownPoints ) const {
         
         std::vector<int> ncomb( 5, 1 );
         std::vector<int> ient( 5, 1 );
         std::vector<T> d( 5, 0.0 );
         
         int kd(0);
         int m(1);
         int ja(0);
         for( int i = 0; i < nDimensions; ++i ) {
            ncomb.at( i ) = 1;
            int jb = ja - 1 + arrayDimensions.at( i );
            
            bool foundJ(false);
            int j(-999);
            for( j = ja; j <= jb; ++j ) {
               if( interpolationPoint.at( i ) <= knownPoints.at( j ) ) {
                  foundJ = true;
                  break;
               } // if
            } // for
            if( not foundJ ) j = jb;
            if( j == ja ) ++j;
            
            d.at( i ) = ( knownPoints.at( j ) - interpolationPoint.at( i ) )
            / ( knownPoints.at( j ) - knownPoints.at( j - 1 ) );
            
            ient.at( i ) = j - ja;
            kd = kd + ient.at( i ) * m;
            m = m * arrayDimensions.at( i );
            ja = jb + 1;
         } // for
         
         T result( 0.0 );
         
         while( true ) {
            T fac = 1.0;
            int iadr = kd;
            int ifadr = 1;
            
            for( int i = 0; i < nDimensions; ++i ) {
               if( ncomb.at( i ) == 0 ) {
                  fac = fac * d.at( i );
                  iadr = iadr - ifadr;
               } // if...
               else {
                  fac = fac * ( 1.0 - d.at( i ) );
               } // else
               ifadr = ifadr * arrayDimensions.at( i );
            } // for i...
            
            result = result + fac * valuesAtKnownPoints.at( iadr );
            int il = nDimensions - 1;
            
            while( ncomb.at( il ) == 0 ) {
               --il;
               if( il == -1 ) {
                  // We found the value at the interpolated point!
                  return result;
               } // if
            } // if
            
            ncomb.at( il ) = 0;
            if( il == nDimensions - 1 ) {
               continue; // Reloop back from top of while(true)
            } // if
            ++il;
            
            for( int k = il; k < nDimensions; ++k ) {
               ncomb.at( k ) = 1;
            } // for( int k...         
            continue; // Reloop back from top of while(true)
         } // if
         
         // I don't think we should ever actually get to this point:
         return result;
      }
      
   };
   
} // namespace tppmc

#endif
